Assumptions
Simulate non-linearity (quadratic)
set.seed(221020)
quadr <- function(x, a, b, c){ a + b*x + c*x^2 }
df <- tibble(
x = seq(-2, 2, length.out = 101),
y_good = abs(4 + (2 * x) + rnorm(101, 0, 1)),
y_quad = abs(quadr(x + 0.5, 1, 0, 1) + rnorm(101, 0, 1)),
)
p_good <- df %>%
ggplot(aes(x=x, y=y_good)) +
geom_point() +
geom_smooth(method = 'lm', se = FALSE) +
NULL
p_quad <- df |>
ggplot(aes(x=x, y=y_quad)) +
geom_point() +
geom_smooth(method = 'lm', se = FALSE) +
NULL
p_good + p_quad

check with component residual plots
component = the individual association of a particular predictor and
the outcome, with all the other predictors held constant.
residual = the differences between the line the model predicts and
the actual observed data points.
df
NA
NA
Simulate autocorrelated errors
If you roll a 6-sided die at 9am everyday, you’ll find there is no
correlation from one day to the next. Every time you roll the die you’ll
get a completely random number coming from the die. If you look at the
temperature at 9am everyday, you will notice that days tend to be
correlated with adjacent days. If it’s 10 degrees on Monday, it’s far
more likely to be 0-20 degrees on Tuesday than it is to be -20 degrees
or 50 degrees. Even though it might be quite common to have a 50 degree
morning, the previous data suggests it probably won’t be.
Probably needlessly complex, esp with so much other stuff to cover.
Worth mentioning conceptually but not belabouring.
# source: https://stackoverflow.com/a/77090029
acf.negbin <- function(N, mu, size, alpha, max.iter = 100, tol = 1e-5) {
m = length(alpha)
generate = function(){
x = sort(rnbinom(N,size=size,mu=mu))
y <- rnorm(length(x))
x[rank(stats::filter(y, alpha, circular = TRUE))]
}
a = generate()
iter <- 0L
ACF <- function(x) acf(x, lag.max = m - 1, plot = FALSE)$acf[1:m]
SSE <- function(x,alpha) sum((ACF(x) - alpha)^2)
while ((SSE0 <- SSE(a, alpha)) > tol) {
if ((iter <- iter + 1L) == max.iter) break
a1 <- generate()
if(SSE(a1,alpha) < SSE0) a <- a1
}
return(a)
}
set.seed(2023)
acf_hi = acf.negbin(
30, # how many obs to simulate
mu = 10, # mean of simulated data
size = 5, # dispersion parameter.
# the bigger, the tighter the data will cluster around the mean.
# alpha=c(1,0.5) # ? correl of points with lag = 1 and lag = 2
alpha=c(0.9, 0.8, 0.8, 0.8) # ? correl of points with lag = 1 and lag = 2
)
acf_hi
[1] 13 5 10 9 5 7 5 4 5 6 7 5 5 6 7 7 10 10 11 15 12 21 17 21 23 16 11 13 9 11
plot(acf_hi,type='l')

# mean(x)
acf(acf_hi)#acf[1:2]

set.seed(2023)
acf_lo <- sample(1:20, size = 30, replace = TRUE)
plot(acf_lo, type='l')

acf(acf_lo)

Simulate non-normal errors
Let’s try generating non-normal errors by sampling from an
asymmetrical beta distribution. Centre those estims on zero and then use
them as errors. So it’ll be normal in that the mean is zero, but they’ll
be super skewed.
set.seed(1)
beta_errors <- rbeta(101, shape1 = 1, shape2 = 50) |>
scale(scale = FALSE, center = TRUE) * 50
plot(beta_errors)

set.seed(221020)
df <- tibble(
x = seq(-2, 2, length.out = 101),
y_good = abs(4 + (2 * x) + rnorm(101, 0, 1)),
y_beta = abs(4 + (2 * x) + beta_errors),
# y_quad = abs(quadr(x + 0.5, 1, 0, 1) + rnorm(101, 0, 1)),
)
# df2
p_good <- df |>
ggplot(aes(x=x, y=y_good)) +
geom_point() +
geom_smooth(method = 'lm', se = FALSE) +
labs(y = 'y') +
NULL
p_beta <- df |>
ggplot(aes(x=x, y=y_beta)) +
geom_point() +
geom_smooth(method = 'lm', se = FALSE) +
labs(y = 'y') +
NULL
p_good + p_beta

Plot normality of errors
m_good <- lm(y_good ~ x, df)
m_beta <- lm(y_beta ~ x, df)
resid_df <- tibble(
good_resid = m_good$residuals,
beta_resid = m_beta$residuals,
)
sd_good <- sd(m_good$residuals)
sd_beta <- sd(m_beta$residuals)
p_good_resid <- resid_df |>
ggplot() +
geom_histogram(aes(x = good_resid)) +
geom_function(
fun = function(x) dnorm(x, mean = 0, sd = sd_good) * 12,
linewidth = 2) +
xlim(-1.5, 1.5) +
labs(
# title = 'normal residuals'
) +
NULL
p_beta_resid <- resid_df |>
ggplot() +
geom_histogram(aes(x = beta_resid)) +
geom_function(
fun = function(x) dnorm(x, mean = 0, sd = .95) * 35,
linewidth = 2
) +
xlim(-4.5, 4.5) +
labs(
# title = 'non-normal residuals'
) +
NULL
p_good_resid + p_beta_resid

q-q plot
plot(m_good, which = 2)

plot(m_beta, which = 2)

Simulate heteroscedasticity
Errors that increase with x
set.seed(221020)
df <- tibble(
x = seq(-2, 2, length.out = 101),
y_good = abs(4 + (2 * x) + rnorm(101, 0, 1)),
y_unequal = abs(4 + (2 * x) +
rnorm(101, mean = 0, sd = seq(0.1, 4, length.out = 101))
),
)
p_unequal <- df |>
ggplot(aes(x=x, y=y_unequal)) +
geom_point() +
geom_smooth(method = 'lm', se = FALSE) +
NULL
p_good + p_unequal

m_unequal <- lm(y_unequal ~ x, data = df)
car::residualPlot(m_unequal) # same as plot() so idk if we need it

plot(m_good, which = 1)

plot(m_unequal, which = 1)

LS0tCnRpdGxlOiAiTGVjdHVyZSAwOCBwbGF5Z3JvdW5kIgpvdXRwdXQ6IAogIGh0bWxfbm90ZWJvb2s6CiAgICB0b2M6IHRydWUKLS0tCgpgYGB7cn0KbGlicmFyeSh0aWR5dmVyc2UpCmxpYnJhcnkoY2FyKQpsaWJyYXJ5KHBlcmZvcm1hbmNlKQpsaWJyYXJ5KHBhdGNod29yaykKYGBgCgoKCiMgQXNzdW1wdGlvbnMKCiMjIFNpbXVsYXRlIG5vbi1saW5lYXJpdHkgKHF1YWRyYXRpYykKCmBgYHtyIG1lc3NhZ2UgPSBGQUxTRX0Kc2V0LnNlZWQoMjIxMDIwKQpxdWFkciA8LSBmdW5jdGlvbih4LCBhLCBiLCBjKXsgYSArIGIqeCArIGMqeF4yICB9CgpkZiA8LSB0aWJibGUoCiAgeCA9IHNlcSgtMiwgMiwgbGVuZ3RoLm91dCA9IDEwMSksCiAgeV9nb29kID0gYWJzKDQgKyAoMiAqIHgpICsgcm5vcm0oMTAxLCAwLCAxKSksCiAgeV9xdWFkID0gYWJzKHF1YWRyKHggKyAwLjUsIDEsIDAsIDEpICsgcm5vcm0oMTAxLCAwLCAxKSksCikKCnBfZ29vZCA8LSBkZiAlPiUKICBnZ3Bsb3QoYWVzKHg9eCwgeT15X2dvb2QpKSArCiAgZ2VvbV9wb2ludCgpICsKICBnZW9tX3Ntb290aChtZXRob2QgPSAnbG0nLCBzZSA9IEZBTFNFKSArCiAgTlVMTAoKcF9xdWFkIDwtIGRmIHw+CiAgZ2dwbG90KGFlcyh4PXgsIHk9eV9xdWFkKSkgKwogIGdlb21fcG9pbnQoKSArCiAgZ2VvbV9zbW9vdGgobWV0aG9kID0gJ2xtJywgc2UgPSBGQUxTRSkgKwogIE5VTEwKCnBfZ29vZCArIHBfcXVhZApgYGAKCiMjIyBjaGVjayB3aXRoIGNvbXBvbmVudCByZXNpZHVhbCBwbG90cwoKY29tcG9uZW50ID0gdGhlIGluZGl2aWR1YWwgYXNzb2NpYXRpb24gb2YgYSBwYXJ0aWN1bGFyIHByZWRpY3RvciBhbmQgdGhlIG91dGNvbWUsIHdpdGggYWxsIHRoZSBvdGhlciBwcmVkaWN0b3JzIGhlbGQgY29uc3RhbnQuCgpyZXNpZHVhbCA9IHRoZSBkaWZmZXJlbmNlcyBiZXR3ZWVuIHRoZSBsaW5lIHRoZSBtb2RlbCBwcmVkaWN0cyBhbmQgdGhlIGFjdHVhbCBvYnNlcnZlZCBkYXRhIHBvaW50cy4KCgpgYGB7cn0KZGYKCgpgYGAKCgoKCiMjIFNpbXVsYXRlIGF1dG9jb3JyZWxhdGVkIGVycm9ycwoKPiBJZiB5b3Ugcm9sbCBhIDYtc2lkZWQgZGllIGF0IDlhbSBldmVyeWRheSwgeW91J2xsIGZpbmQgdGhlcmUgaXMgbm8gY29ycmVsYXRpb24gZnJvbSBvbmUgZGF5IHRvIHRoZSBuZXh0LiBFdmVyeSB0aW1lIHlvdSByb2xsIHRoZSBkaWUgeW91J2xsIGdldCBhIGNvbXBsZXRlbHkgcmFuZG9tIG51bWJlciBjb21pbmcgZnJvbSB0aGUgZGllLgpJZiB5b3UgbG9vayBhdCB0aGUgdGVtcGVyYXR1cmUgYXQgOWFtIGV2ZXJ5ZGF5LCB5b3Ugd2lsbCBub3RpY2UgdGhhdCBkYXlzIHRlbmQgdG8gYmUgY29ycmVsYXRlZCB3aXRoIGFkamFjZW50IGRheXMuIElmIGl0J3MgMTAgZGVncmVlcyBvbiBNb25kYXksIGl0J3MgZmFyIG1vcmUgbGlrZWx5IHRvIGJlIDAtMjAgZGVncmVlcyBvbiBUdWVzZGF5IHRoYW4gaXQgaXMgdG8gYmUgLTIwIGRlZ3JlZXMgb3IgNTAgZGVncmVlcy4gRXZlbiB0aG91Z2ggaXQgbWlnaHQgYmUgcXVpdGUgY29tbW9uIHRvIGhhdmUgYSA1MCBkZWdyZWUgbW9ybmluZywgdGhlIHByZXZpb3VzIGRhdGEgc3VnZ2VzdHMgaXQgcHJvYmFibHkgd29uJ3QgYmUuCgpQcm9iYWJseSBuZWVkbGVzc2x5IGNvbXBsZXgsIGVzcCB3aXRoIHNvIG11Y2ggb3RoZXIgc3R1ZmYgdG8gY292ZXIuCldvcnRoIG1lbnRpb25pbmcgY29uY2VwdHVhbGx5IGJ1dCBub3QgYmVsYWJvdXJpbmcuCgoKYGBge3J9CiMgc291cmNlOiBodHRwczovL3N0YWNrb3ZlcmZsb3cuY29tL2EvNzcwOTAwMjkgCgphY2YubmVnYmluIDwtIGZ1bmN0aW9uKE4sIG11LCBzaXplLCBhbHBoYSwgbWF4Lml0ZXIgPSAxMDAsIHRvbCA9IDFlLTUpIHsKICAKICBtID0gbGVuZ3RoKGFscGhhKQogIAogIGdlbmVyYXRlID0gZnVuY3Rpb24oKXsKICAgeCA9IHNvcnQocm5iaW5vbShOLHNpemU9c2l6ZSxtdT1tdSkpCiAgIHkgPC0gcm5vcm0obGVuZ3RoKHgpKQogICB4W3Jhbmsoc3RhdHM6OmZpbHRlcih5LCBhbHBoYSwgY2lyY3VsYXIgPSBUUlVFKSldCiAgfQogIAogIGEgPSBnZW5lcmF0ZSgpCiAgaXRlciA8LSAwTAogIAogIEFDRiA8LSBmdW5jdGlvbih4KSBhY2YoeCwgbGFnLm1heCA9IG0gLSAxLCBwbG90ID0gRkFMU0UpJGFjZlsxOm1dCiAgU1NFIDwtIGZ1bmN0aW9uKHgsYWxwaGEpIHN1bSgoQUNGKHgpIC0gYWxwaGEpXjIpCiAgICAKICB3aGlsZSAoKFNTRTAgPC0gU1NFKGEsIGFscGhhKSkgPiB0b2wpIHsKICAgIGlmICgoaXRlciA8LSBpdGVyICsgMUwpID09IG1heC5pdGVyKSBicmVhawogICAgYTEgPC0gZ2VuZXJhdGUoKQogICAgaWYoU1NFKGExLGFscGhhKSA8IFNTRTApIGEgPC0gYTEKICB9CiAgICAKICByZXR1cm4oYSkKCn0KYGBgCgpgYGB7cn0Kc2V0LnNlZWQoMjAyMykKYWNmX2hpID0gYWNmLm5lZ2JpbigKICAzMCwgICAgICAgICAgICAjIGhvdyBtYW55IG9icyB0byBzaW11bGF0ZQogIG11ID0gMTAsICAgICAgIyBtZWFuIG9mIHNpbXVsYXRlZCBkYXRhCiAgc2l6ZSA9IDUsICAgICAjIGRpc3BlcnNpb24gcGFyYW1ldGVyLiAKICAgICAgICAgICAgICAgICAgIyB0aGUgYmlnZ2VyLCB0aGUgdGlnaHRlciB0aGUgZGF0YSB3aWxsIGNsdXN0ZXIgYXJvdW5kIHRoZSBtZWFuLgogICMgYWxwaGE9YygxLDAuNSkgICMgPyBjb3JyZWwgb2YgcG9pbnRzIHdpdGggbGFnID0gMSBhbmQgbGFnID0gMiAKICBhbHBoYT1jKDAuOSwgMC44LCAwLjgsIDAuOCkgICMgPyBjb3JyZWwgb2YgcG9pbnRzIHdpdGggbGFnID0gMSBhbmQgbGFnID0gMiAKICApCgphY2ZfaGkKYGBgCgoKYGBge3J9CnBsb3QoYWNmX2hpLHR5cGU9J2wnKQojIG1lYW4oeCkKYWNmKGFjZl9oaSkjYWNmWzE6Ml0KYGBgCgpgYGB7cn0Kc2V0LnNlZWQoMjAyMykKYWNmX2xvIDwtIHNhbXBsZSgxOjIwLCBzaXplID0gMzAsIHJlcGxhY2UgPSBUUlVFKQoKcGxvdChhY2ZfbG8sIHR5cGU9J2wnKQphY2YoYWNmX2xvKQpgYGAKCgoKIyMgU2ltdWxhdGUgbm9uLW5vcm1hbCBlcnJvcnMKCkxldCdzIHRyeSBnZW5lcmF0aW5nIG5vbi1ub3JtYWwgZXJyb3JzIGJ5IHNhbXBsaW5nIGZyb20gYW4gYXN5bW1ldHJpY2FsIGJldGEgZGlzdHJpYnV0aW9uLgpDZW50cmUgdGhvc2UgZXN0aW1zIG9uIHplcm8gYW5kIHRoZW4gdXNlIHRoZW0gYXMgZXJyb3JzLgpTbyBpdCdsbCBiZSBub3JtYWwgaW4gdGhhdCB0aGUgbWVhbiBpcyB6ZXJvLCBidXQgdGhleSdsbCBiZSBzdXBlciBza2V3ZWQuCgpgYGB7cn0Kc2V0LnNlZWQoMSkKYmV0YV9lcnJvcnMgPC0gcmJldGEoMTAxLCBzaGFwZTEgPSAxLCBzaGFwZTIgPSA1MCkgfD4gCiAgc2NhbGUoc2NhbGUgPSBGQUxTRSwgY2VudGVyID0gVFJVRSkgKiA1MApwbG90KGJldGFfZXJyb3JzKQpgYGAKCgpgYGB7ciBtZXNzYWdlID0gRn0Kc2V0LnNlZWQoMjIxMDIwKQoKZGYgPC0gdGliYmxlKAogIHggPSBzZXEoLTIsIDIsIGxlbmd0aC5vdXQgPSAxMDEpLAogIHlfZ29vZCA9IGFicyg0ICsgKDIgKiB4KSArIHJub3JtKDEwMSwgMCwgMSkpLCAKICB5X2JldGEgPSBhYnMoNCArICgyICogeCkgKyBiZXRhX2Vycm9ycyksCiAgCiAgIyB5X3F1YWQgPSBhYnMocXVhZHIoeCArIDAuNSwgMSwgMCwgMSkgKyBybm9ybSgxMDEsIDAsIDEpKSwKKQojIGRmMgoKcF9nb29kIDwtIGRmIHw+CiAgZ2dwbG90KGFlcyh4PXgsIHk9eV9nb29kKSkgKwogIGdlb21fcG9pbnQoKSArCiAgZ2VvbV9zbW9vdGgobWV0aG9kID0gJ2xtJywgc2UgPSBGQUxTRSkgKwogIGxhYnMoeSA9ICd5JykgKwogIE5VTEwKCnBfYmV0YSA8LSBkZiB8PgogIGdncGxvdChhZXMoeD14LCB5PXlfYmV0YSkpICsKICBnZW9tX3BvaW50KCkgKwogIGdlb21fc21vb3RoKG1ldGhvZCA9ICdsbScsIHNlID0gRkFMU0UpICsKICBsYWJzKHkgPSAneScpICsKICBOVUxMCgpwX2dvb2QgKyBwX2JldGEKYGBgCgoKIyMjIFBsb3Qgbm9ybWFsaXR5IG9mIGVycm9ycwoKYGBge3J9Cm1fZ29vZCA8LSBsbSh5X2dvb2QgfiB4LCBkZikKbV9iZXRhIDwtIGxtKHlfYmV0YSB+IHgsIGRmKQpgYGAKCmBgYHtyIG1lc3NhZ2UgPSBGfQpyZXNpZF9kZiA8LSB0aWJibGUoCiAgZ29vZF9yZXNpZCA9IG1fZ29vZCRyZXNpZHVhbHMsCiAgYmV0YV9yZXNpZCA9IG1fYmV0YSRyZXNpZHVhbHMsCikKCnNkX2dvb2QgPC0gc2QobV9nb29kJHJlc2lkdWFscykKc2RfYmV0YSA8LSBzZChtX2JldGEkcmVzaWR1YWxzKQoKcF9nb29kX3Jlc2lkIDwtIHJlc2lkX2RmIHw+CiAgZ2dwbG90KCkgKwogIGdlb21faGlzdG9ncmFtKGFlcyh4ID0gZ29vZF9yZXNpZCkpICsKICBnZW9tX2Z1bmN0aW9uKAogICAgZnVuID0gZnVuY3Rpb24oeCkgZG5vcm0oeCwgbWVhbiA9IDAsIHNkID0gc2RfZ29vZCkgKiAxMiwKICAgIGxpbmV3aWR0aCA9IDIpICsKICB4bGltKC0xLjUsIDEuNSkgKwogIGxhYnMoCiAgICAjIHRpdGxlID0gJ25vcm1hbCByZXNpZHVhbHMnCiAgKSArCiAgTlVMTAoKcF9iZXRhX3Jlc2lkIDwtIHJlc2lkX2RmIHw+CiAgZ2dwbG90KCkgKwogIGdlb21faGlzdG9ncmFtKGFlcyh4ID0gYmV0YV9yZXNpZCkpICsKICBnZW9tX2Z1bmN0aW9uKAogICAgZnVuID0gZnVuY3Rpb24oeCkgZG5vcm0oeCwgbWVhbiA9IDAsIHNkID0gLjk1KSAqIDM1LAogICAgbGluZXdpZHRoID0gMgogICkgKwogIHhsaW0oLTQuNSwgNC41KSArCiAgbGFicygKICAgICMgdGl0bGUgPSAnbm9uLW5vcm1hbCByZXNpZHVhbHMnCiAgKSArCiAgTlVMTAoKcF9nb29kX3Jlc2lkICsgcF9iZXRhX3Jlc2lkCmBgYAoKIyMjIHEtcSBwbG90CgpgYGB7cn0KcGxvdChtX2dvb2QsIHdoaWNoID0gMikKcGxvdChtX2JldGEsIHdoaWNoID0gMikKYGBgCgoKCgojIyBTaW11bGF0ZSBoZXRlcm9zY2VkYXN0aWNpdHkKCkVycm9ycyB0aGF0IGluY3JlYXNlIHdpdGggeAoKYGBge3IgbWVzc2FnZSA9IEZ9CnNldC5zZWVkKDIyMTAyMCkKCmRmIDwtIHRpYmJsZSgKICB4ID0gc2VxKC0yLCAyLCBsZW5ndGgub3V0ID0gMTAxKSwKICB5X2dvb2QgPSBhYnMoNCArICgyICogeCkgKyBybm9ybSgxMDEsIDAsIDEpKSwKICB5X3VuZXF1YWwgPSBhYnMoNCArICgyICogeCkgKyAKICAgICAgICAgICAgICAgICAgICBybm9ybSgxMDEsIG1lYW4gPSAwLCBzZCA9IHNlcSgwLjEsIDQsIGxlbmd0aC5vdXQgPSAxMDEpKQogICAgICAgICAgICAgICAgICApLAopCgpwX3VuZXF1YWwgPC0gZGYgfD4KICBnZ3Bsb3QoYWVzKHg9eCwgeT15X3VuZXF1YWwpKSArCiAgZ2VvbV9wb2ludCgpICsKICBnZW9tX3Ntb290aChtZXRob2QgPSAnbG0nLCBzZSA9IEZBTFNFKSArCiAgTlVMTAoKcF9nb29kICsgcF91bmVxdWFsCmBgYAoKYGBge3J9Cm1fdW5lcXVhbCA8LSBsbSh5X3VuZXF1YWwgfiB4LCBkYXRhID0gZGYpCgpjYXI6OnJlc2lkdWFsUGxvdChtX3VuZXF1YWwpICAjIHNhbWUgYXMgcGxvdCgpIHNvIGlkayBpZiB3ZSBuZWVkIGl0CmBgYAoKYGBge3J9CnBsb3QobV9nb29kLCB3aGljaCA9IDEpCnBsb3QobV91bmVxdWFsLCB3aGljaCA9IDEpCmBgYAoKCiMgRGlhZ25vc3RpY3MKCiMjIFRocmVlIGNoZWNrcyBmb3IgdW51c3VhbCBkYXRhIHBvaW50cwoKTWFrZSBzdXJlIHNpbmdsZSBwb2ludHMgYXJlbid0IGRyaXZpbmcgdGhlIHJlc3VsdHMuCgojIyMgT3V0bGllcnMKCkV4dHJlbWUgdmFsdWVzIG9mIHkuCgpIb3cgcXVhbnRpZnk/ClVuZXhwZWN0ZWRseSBsYXJnZSAoc3R1ZGVudGlzZWQpIHJlc2lkdWFscy4KCiMjIyMgU2ltdWxhdGUKCmBgYHtyIG1lc3NhZ2UgPSBGfQpzZXQuc2VlZCgyMjEwMjApCmRmIDwtIHRpYmJsZSgKICB4ID0gc2VxKC0yLCAyLCBsZW5ndGgub3V0ID0gMTAxKSwKICB5X2dvb2QgPSBhYnMoNCArICgyICogeCkgKyBybm9ybSgxMDEsIDAsIDEpKSwKICB5X291dGwgPSB5X2dvb2QKKQoKZGZbOTUsICd5X291dGwnXSA8LSAxLjIzNjIyNwojIGRmWzk4LCAneV9vdXRsJ10gPC0gMi42MjM5NzMKIyBkZlszMCwgJ3lfb3V0bCddIDwtIDUuNjQzMDMzCiMgZGZbMjAsICd5X291dGwnXSA8LSAzLjE2MDI2MQojIGRmWzE0LCAneV9vdXRsJ10gPC0gMy43NjIzNTIKCnBfb3V0bCA8LSBkZiB8PgogIGdncGxvdChhZXMoeD14LCB5PXlfb3V0bCkpICsKICBnZW9tX3BvaW50KCkgKwogIGdlb21fc21vb3RoKG1ldGhvZCA9ICdsbScsIHNlID0gRkFMU0UpICsKICBOVUxMCgpwX2dvb2QgKyBwX291dGwKYGBgCgoKCgojIyMjIENoZWNrIHdpdGggc3R1ZGVudGlzZWQgcmVzaWR1YWxzCgp1bnN0YW5kYXJkaXNlZCByZXNpZHVhbHMgLT4gIHN0YW5kYXJkaXNlZCByZXNpZHVhbHMgIC0+IHN0dWRlbnRpc2VkIHJlc2lkdWFscyB3aXRoIGByc3R1ZGVudChtb2RlbClgCgpgYGB7cn0KbV9vdXRsIDwtIGxtKHlfb3V0bCB+IHgsIGRhdGEgPSBkZikKCmRmJHN0dWRfcmVzaWQgPC0gcnN0dWRlbnQobV9vdXRsKQoKZGYgfD4KICBnZ3Bsb3QoYWVzKHggPSB4LCB5ID0gc3R1ZF9yZXNpZCkpICsKICBnZW9tX3BvaW50KCkgKwogIGdlb21faGxpbmUoeWludGVyY2VwdCA9ICAyLCBjb2xvdXIgPSAncmVkJykgKwogIGdlb21faGxpbmUoeWludGVyY2VwdCA9IC0yLCBjb2xvdXIgPSAncmVkJykgKwogIE5VTEwKCmRmIHw+CiAgZmlsdGVyKHN0dWRfcmVzaWQgPiAyIHwgc3R1ZF9yZXNpZCA8IC0yKQpgYGAKCgoKIyMjIEhpZ2ggbGV2ZXJhZ2UKCiMjIyMgU2ltdWxhdGUKCmBgYHtyIG1lc3NhZ2UgPSBGfQpzZXQuc2VlZCgyMjEwMjApCgpkZiA8LSB0aWJibGUoCiAgeCA9IHNlcSgtMiwgMiwgbGVuZ3RoLm91dCA9IDEwMyksCiAgeV9nb29kID0gYWJzKDQgKyAoMiAqIHgpICsgcm5vcm0oMTAzLCAwLCAxKSksCiAgeV9sZXZnID0gYWJzKDQgKyAoMiAqIHgpICsgcm5vcm0oMTAzLCAwLCAxKSksCikKCgpkZiA8LSByYmluZCgKICBkZiwKICB0aWJibGUoCiAgICB4ID0gYygKICAgICAgMwogICAgICAgICAgIyA0LjUKICAgICAgICAgICMgLTMuNQogICAgICAgICAgKSwKICAgIHlfZ29vZCA9IE5BLAogICAgeV9sZXZnID0gYygKICAgICAgNi4wNjgxNzIKICAgICAgIyAxMi4yMDE4MjcKICAgICAgIyAyLjIzNjgwMzAKICAgICAgKQogICkKKQogIApwX2xldmcgPC0gZGYgfD4KICBnZ3Bsb3QoYWVzKHg9eCwgeT15X2xldmcpKSArCiAgZ2VvbV9wb2ludCgpICsKICBnZW9tX3Ntb290aChtZXRob2QgPSAnbG0nLCBzZSA9IEZBTFNFKSArCiAgTlVMTAoKIyBwX2dvb2QgKyBwX2xldmcKcF9sZXZnCmBgYAoKCiMjIyMgQ2hlY2sgd2l0aCBoYXQgdmFsdWVzCgpoYXQgcmVmZXJzIHRvIHRoZSBjaXJjdW1mbGV4IG9uIHRvcCBvZiBhIHZhcmlhYmxlLCB3aGljaCBpcyBsaWtlIGEgbGl0dGxlIHdheSBvZiBzYXlpbmcgInRoaXMgaXMgYW4gZXN0aW1hdGUgZnJvbSBhIG1vZGVsIi4KaXQncyBsaWtlIGluIG5hdHVyYWwgbGFuZ3VhZ2UsIGxpa2UgaW4gRnJlbmNoIGhvdyBhIGNpcmN1bWZsZXggbWVhbnMgInRoZXJlIHVzZWQgdG8gYmUgYW4gJ3MnIGhlcmUiLgpqdXN0IGFuIGV4dHJhIGFubm90YXRpb24uCgpgaGF0dmFsdWVzKClgIGlzIG91ciB3b3JraG9yc2UgaGVyZS4KCmBgYHtyfQojIGZpdCBtb2RlbAptX2xldmcgPC0gbG0oeV9sZXZnIH4geCwgZGF0YSA9IGRmKQoKIyBmaW5kIG1lYW4gZXhwZWN0ZWQgaGF0IHZhbHVlIApuX3ByZWRpY3RvcnMgPC0gMQpuX29ic2VydmF0aW9ucyA8LSBucm93KGRmKQptZWFuX2hhdCA8LSAobl9wcmVkaWN0b3JzKzEpL25fb2JzZXJ2YXRpb25zCgojIHBvaW50cyA+IDIqbWVhbl9oYXQgd2lsbCBiZSBjb25zaWRlcmVkIGhpZ2gtbGV2ZXJhZ2UKZGYkaGF0dmFscyA8LSBoYXR2YWx1ZXMobV9sZXZnKQpkZiB8PgogIGdncGxvdChhZXMoeCA9IHgsIHkgPSBoYXR2YWxzKSkgKwogIGdlb21fcG9pbnQoKSArCiAgZ2VvbV9obGluZSh5aW50ZXJjZXB0ID0gMiptZWFuX2hhdCwgY29sb3VyID0gJ3JlZCcpCgojIFBsb3QgZ29vZCBoYXQgdmFsdWVzCnRpYmJsZSgKICB4ID0gc2VxKC0yLCAyLCBsZW5ndGgub3V0ID0gMTAxKSwKICBoYXR2YWxzID0gaGF0dmFsdWVzKG1fZ29vZCkKKSB8PgogIGdncGxvdChhZXMoeCA9IHgsIHkgPSBoYXR2YWxzKSkgKwogIGdlb21fcG9pbnQoKSArCiAgZ2VvbV9obGluZSh5aW50ZXJjZXB0ID0gMiooMi8xMDEpLCBjb2xvdXIgPSAncmVkJykKYGBgCgpUaGUgY3VydmUgc21pbGV5IGZhY2UgaXMgYmVjYXVzZSB3ZSdyZSBsb29raW5nIGF0IGhvdyBmYXIgdGhlIHZhbHVlcyBhcmUgZnJvbSB0aGUgbWVhbiBvZiB0aGUgZGF0YSB3aGVuIHggPSAwPwpUaGUgZmFydGhlciB5b3UgZ2V0IGZyb20gdGhlIHkgdmFsdWUgd2hlbiB4PTAsIHRoZSBtb3JlIHRoZSBoYXQgdmFsdWVzIHdpbGwgaW5jcmVhc2U/CgpUaGUgcGxvdCBvZiB0aGUgZ29vZCBkYXRhIGFsc28gZ29lcyB0byBzaG93IHRoYXQgd2Ugc2hvdWxkbid0IGF1dG9tYXRpY2FsbHkgYmUgc3VzcGljaW91cy4KSnVzdCBpZiBpdCdzIHJlYWxseSBvdXRzaWRlIG9mIHRoZSBjdXJ2ZSB0aGF0IHdlJ2QgZXhwZWN0LgpEaXNjcmV0aW9uISAKTm90IGp1c3QgYmluYXJ5IGRlY2lzaW9ucy4KCgojIyMgSGlnaCBpbmZsdWVuY2UKCiMjIyMgU2ltdWxhdGUKCmBgYHtyIG1lc3NhZ2UgPSBGfQpzZXQuc2VlZCgyMjEwMjApCgpkZiA8LSB0aWJibGUoCiAgeCA9IHNlcSgtMiwgMiwgbGVuZ3RoLm91dCA9IDEwMSksCiAgeV9nb29kID0gYWJzKDQgKyAoMiAqIHgpICsgcm5vcm0oMTAxLCAwLCAxKSksCiAgeV9pbmZsID0gYWJzKDQgKyAoMiAqIHgpICsgcm5vcm0oMTAxLCAwLCAxKSksCikKCiMgVXNlIHRoZSBzYW1lIG91dGxpZXIgdmFsdWVzCmRmWzk1LCAneV9pbmZsJ10gPC0gMS4yMzYyMjcKZGZbOTgsICd5X2luZmwnXSA8LSAyLjYyMzk3MwpkZlszMCwgJ3lfaW5mbCddIDwtIDUuNjQzMDMzCmRmWzIwLCAneV9pbmZsJ10gPC0gMy4xNjAyNjEKZGZbMTQsICd5X2luZmwnXSA8LSAzLjc2MjM1MgoKIyBBZGQgaW4gdGhlIHNhbWUgaGlnaC1sZXZlcmFnZSB2YWx1ZXMKZGYgPC0gcmJpbmQoCiAgZGYsCiAgdGliYmxlKAogICAgeCA9IGMoMywKICAgICAgICAgIDQuNSwKICAgICAgICAgIC0zLjUpLAogICAgeV9nb29kID0gTkEsCiAgICB5X2luZmwgPSBjKAogICAgICA2LjA2ODE3MiwKICAgICAgMTIuMjAxODI3LAogICAgICAyLjIzNjgwMzApCiAgKQopCgpwX2luZmwgPC0gZGYgfD4KICBnZ3Bsb3QoYWVzKHg9eCwgeT15X2luZmwpKSArCiAgZ2VvbV9wb2ludCgpICsKICBnZW9tX3Ntb290aChtZXRob2QgPSAnbG0nLCBzZSA9IEZBTFNFKSArCiAgTlVMTAoKcF9nb29kICsgcF9pbmZsCmBgYAoKCiMjIyMgQ2hlY2sgd2l0aCBDb29rJ3MgZGlzdGFuY2UKCkNvb2sncyBkaXN0YW5jZSBpcyBiYXNpY2FsbHkgb3V0bHlpbmduZXNzICogbGV2ZXJhZ2UuClNvIGl0J3MgYmlnIGlmIGVpdGhlciBvbmUgb2YgdGhvc2UgYXJlIHRydWUsIGFuZCB2ZXJ5IGJpZyBpZiBib3RoIGFyZSB0cnVlLgoKSW50ZXJwcmV0YXRpb246IHRoZSBhdmVyYWdlIGRpc3RhbmNlIHRoZSBwcmVkaWN0ZWQgb3V0Y29tZSB2YWx1ZXMgd2lsbCBtb3ZlLCBpZiBhIGdpdmVuIGNhc2UgaXMgcmVtb3ZlZC4KClBvc3NpYmxlIHRocmVzaG9sZHM6CgotIDEKLSA0IC8gKHNhbXBsZV9zaXplIC0gbl9wYXJhbXMgLSAxKQotIHJlbGF0aXZlIHRvIGFsbCBEcyBpbiB0aGUgZGF0YSBzZXQKCmBgYHtyfQptX2luZmwgPC0gbG0oeV9pbmZsIH4geCwgZGF0YSA9IGRmKQpkZiREIDwtIGNvb2tzLmRpc3RhbmNlKG1faW5mbCkKCmRmIHw+CiAgZ2dwbG90KGFlcyh4ID0geCwgeSA9IEQpKSArCiAgZ2VvbV9wb2ludCgpICsKICBnZW9tX2hsaW5lKHlpbnRlcmNlcHQgPSA0IC8gKG5fb2JzZXJ2YXRpb25zIC0gbl9wcmVkaWN0b3JzIC0gMSksIGNvbG91ciA9ICdyZWQnKSArCiAgTlVMTAoKdGliYmxlKAogIHggPSBzZXEoLTIsIDIsIGxlbmd0aC5vdXQgPSAxMDEpLAogIEQgPSBjb29rcy5kaXN0YW5jZShtX2dvb2QpCikgfD4KICBnZ3Bsb3QoYWVzKHggPSB4LCB5ID0gRCkpICsKICBnZW9tX3BvaW50KCkgKwogIGdlb21faGxpbmUoeWludGVyY2VwdCA9IDQgLyAoMTAxIC0gbl9wcmVkaWN0b3JzIC0gMSksIGNvbG91ciA9ICdyZWQnKQpgYGAKCmBgYHtyfQpwbG90KG1faW5mbCwgd2hpY2ggPSA0KQpgYGAKCgoKQ2hlY2sgdGhlIG91dGxpZXIgYW5kIGxldmVyYWdlIGRhdGFzZXRzIHRvbzoKCmBgYHtyfQpwbG90KG1fb3V0bCwgd2hpY2ggPSA0KQpwbG90KG1fbGV2Zywgd2hpY2ggPSA0KQpgYGAKCiMjIyMgQ2hlY2sgd2l0aCBjb3ZhcmlhbmNlIHJhdGlvCgpDb3ZhcmlhbmNlIHJhdGlvIGluZGljYXRlcyB0aGUgaW5mbHVlbmNlIG9mIGEgZGF0YSBwb2ludCBvbiB0aGUgKGNvKXZhcmlhbmNlcyBvZiB0aGUgcmVncmVzc2lvbiBwYXJhbWV0ZXJzLgoKYGBge3J9CmRmJGNvdi5yIDwtIGNvdnJhdGlvKG1faW5mbCkKCmNvdnJfdXBwZXIgPC0gMSArICgoMyAqIChuX3ByZWRpY3RvcnMgKyAxKSkvbl9vYnNlcnZhdGlvbnMpCmNvdnJfbG93ZXIgPC0gMSAtICgoMyAqIChuX3ByZWRpY3RvcnMgKyAxKSkvbl9vYnNlcnZhdGlvbnMpCgpkZiB8PgogIGdncGxvdChhZXMoeCA9IHgsIHkgPSBjb3YucikpICsKICBnZW9tX3BvaW50KCkgKwogIGdlb21faGxpbmUoeWludGVyY2VwdCA9IGNvdnJfdXBwZXIsIGNvbG91ciA9ICdyZWQnKSArCiAgZ2VvbV9obGluZSh5aW50ZXJjZXB0ID0gY292cl9sb3dlciwgY29sb3VyID0gJ3JlZCcpICsKICBOVUxMCgp0aWJibGUoCiAgeCA9IHNlcSgtMiwgMiwgbGVuZ3RoLm91dCA9IDEwMSksCiAgY292LnIgPSBjb3ZyYXRpbyhtX2dvb2QpCikgfD4KICBnZ3Bsb3QoYWVzKHggPSB4LCB5ID0gY292LnIpKSArCiAgZ2VvbV9wb2ludCgpICsKICBnZW9tX2hsaW5lKHlpbnRlcmNlcHQgPSBjb3ZyX3VwcGVyLCBjb2xvdXIgPSAncmVkJykgKwogIGdlb21faGxpbmUoeWludGVyY2VwdCA9IGNvdnJfbG93ZXIsIGNvbG91ciA9ICdyZWQnKSArCiAgTlVMTApgYGAKCkl0IHdpbGwgZGV0ZWN0IHdlaXJkbHkgaGlnaC1pbmZsdWVuY2UgcG9pbnRzIGlmIHRoZXkncmUgdGhlcmUsIGJ1dCBpdCB3aWxsIGFsc28gZGV0ZWN0IGhpZ2gtaW5mbHVlbmNlIHBvaW50cyByZWdhcmRsZXNzLgpKdXN0IGJlY2F1c2Ugc29tZXRoaW5nIGV4Y2VlZHMgdGhlc2UgdGhyZXNob2xkcyBET0VTIE5PVCBNRUFOIElUJ1MgV0VJUkQgT1IgQ09NRVMgRlJPTSBBIERJRkYgRElTVFJJQlVUSU9OL0dFTkVSQVRJVkUgUFJPQ0VTUy4KVGhpcyBpcyB3aHkgd2UgbmVlZCBkaXNjcmV0aW9uLgoKCiMjIyBDaGVjayBhbGwgdGhlc2UgZGlhZ25vc3RpYyBtZWFzdXJlcyBhdCBvbmNlIHdpdGggdGhlIGZ1bmN0aW9uIHRoYXQgZG9lcyBpdCBhbGwKCmBgYHtyfQppbmZsdWVuY2UubWVhc3VyZXMobV9pbmZsKQpgYGAKCi0gYGRmYi4xX2A6IGRpZmZlcmVuY2UgYmV0d2VlbiB0aGUgcHJlZGljdGVkIHZhbHVlcyBmb3IgdGhlIGludGVyY2VwdCB3aXRoIGFuZCB3aXRob3V0IHRoaXMgb2JzZXJ2YXRpb24KLSBgZGZiLnhgOiBkaWZmZXJlbmNlIGJldHdlZW4gdGhlIHByZWRpY3RlZCB2YWx1ZXMgZm9yIHRoZSBzbG9wZSB3aXRoIGFuZCB3aXRob3V0IHRoaXMgb2JzZXJ2YXRpb24gPC0gdGhlcmUgd2lsbCBiZSBvbmUgb2YgdGhlc2UgcGVyIHByZWRpY3RvcgotIGBkZmZpdGA6IGRpZmZlcmVuY2UgYmV0d2VlbiBwcmVkaWN0ZWQgdmFsdWVzIGZvciB0aGUgb3V0Y29tZSB3aXRoIGFuZCB3aXRob3V0IHRoaXMgb2JzZXJ2YXRpb24KLSBgY292LnJgOiByYXRpbyBvZiB0aGUgY292YXJpYW5jZSBvZiByZWdyZXNzaW9uIHBhcmFtZXRlcnMgd2l0aCBhbmQgd2l0aG91dCB0aGlzIG9ic2VydmF0aW9uCi0gYGNvb2suZGA6IENvb2sncyBEaXN0YW5jZSBvZiB0aGlzIG9ic2VydmF0aW9uCi0gYGhhdGA6IGhhdCB2YWx1ZSBvZiB0aGlzIG9ic2VydmF0aW9uCgoKCgoKIyMgRmluYWxseTogb25lIGNoZWNrIGZvciByZWxhdGlvbnNoaXBzIGJldHdlZW4gcHJlZGljdG9ycwoKIyMjIE11bHRpY29sbGluZWFyaXR5CgojIyMjIFNpbXVsYXRlCgpgYGB7cn0Kc2V0LnNlZWQoMSkKCm5fb2JzIDwtIDEwMQpjb3JyX3gxeSAgPC0gMC44CmNvcnJfeDJ5ICA8LSAwLjYKIyBjb3JyX3gxeDIgPC0gMC45ICAjIHZpZiA9IDUuMwpjb3JyX3gxeDIgPC0gMC45MyAgICMgdmlmID0gVEJECgpjb3JyX2RmIDwtIE1BU1M6Om12cm5vcm0oCiAgbiA9IG5fb2JzLAogIG11ID0gYygwLCAwLCAwKSwKICBTaWdtYSA9IG1hdHJpeCgKICAgIGMoMSwgY29ycl94MXksIGNvcnJfeDJ5LAogICAgIGNvcnJfeDF5LCAxLCBjb3JyX3gxeDIsCiAgICAgY29ycl94MnksIGNvcnJfeDF4MiwgMSksCiAgICBucm93ID0gMyksCiAgZW1waXJpY2FsID0gVFJVRQogICkKCmNvbG5hbWVzKGNvcnJfZGYpIDwtIGMoJ3knLCAneDEnLCAneDInKQpjb3JyX2RmIDwtIGFzX3RpYmJsZShjb3JyX2RmKQoKcGFpcnMoY29ycl9kZikKCmNvcnJfZGYgfD4gYXNfdGliYmxlKCkgfD4KICBnZ3Bsb3QoYWVzKHggPSB4MSwgeSA9IHgyKSkgKwogIGdlb21fcG9pbnQoKQpgYGAKCmBgYHtyfQpzZXQuc2VlZCgxKQoKdW5jb3JyX3gxeDIgPC0gMC4wICAjIGlkZWFsbHkgdGhpcyBxdWFudGl0eSB3b3VsZCBiZSBsb3cKCnVuY29ycl9kZiA8LSBNQVNTOjptdnJub3JtKAogIG4gPSBuX29icywKICBtdSA9IGMoMCwgMCwgMCksCiAgU2lnbWEgPSBtYXRyaXgoCiAgICBjKDEsIGNvcnJfeDF5LCBjb3JyX3gyeSwKICAgICBjb3JyX3gxeSwgMSwgdW5jb3JyX3gxeDIsCiAgICAgY29ycl94MnksIHVuY29ycl94MXgyLCAxKSwKICAgIG5yb3cgPSAzKSwKICBlbXBpcmljYWwgPSBUUlVFCiAgKQoKY29sbmFtZXModW5jb3JyX2RmKSA8LSBjKCd5JywgJ3gxJywgJ3gyJykKdW5jb3JyX2RmIDwtIGFzX3RpYmJsZSh1bmNvcnJfZGYpCgpwYWlycyh1bmNvcnJfZGYpCgp1bmNvcnJfZGYgfD4gCiAgZ2dwbG90KGFlcyh4ID0geDEsIHkgPSB4MikpICsKICBnZW9tX3BvaW50KCkKYGBgCgpgYGB7cn0KY29ycl9kZiB8PiAKICBnZ3Bsb3QoYWVzKHggPSB4MSwgeSA9IHksIGNvbG91ciA9IHgyKSkgKwogIGdlb21fcG9pbnQoKQoKdW5jb3JyX2RmIHw+IAogIGdncGxvdChhZXMoeCA9IHgxLCB5ID0geSwgY29sb3VyID0geDIpKSArCiAgZ2VvbV9wb2ludCgpCmBgYAoKCiMjIyMgRml0IG1vZGVscwoKYGBge3J9Cm1fY29yciA8LSBsbSh5IH4geDEgKyB4MiwgZGF0YSA9IGNvcnJfZGYpCnN1bW1hcnkobV9jb3JyKQpgYGAKCgpgYGB7cn0KbV91bmNvcnIgPC0gbG0oeSB+IHgxICsgeDIsIGRhdGEgPSB1bmNvcnJfZGYpCnN1bW1hcnkobV91bmNvcnIpCmBgYAoKCgojIyMjIENoZWNrIHdpdGggVklGCgpgYGB7cn0KdmlmKG1fdW5jb3JyKQpgYGAKCgpgYGB7cn0KdmlmKG1fY29ycikKYGBgCgoKPDUgaXMgbG93LgpiZXR3ZWVuIDUgYW5kIDEwIGlzIG1vZGVyYXRlLgpNb3JlIHRoYW4gMTAgaXMgYSBiaWcgcHJvYmxlbS4KCgoKCiMjIFBsb3RzIHdpdGggYHBlcmZvcm1hbmNlOjpjaGVja19tb2RlbCgpYAoKU29tZSBvZiB0aGVzZSBhc3N1bXB0aW9ucy9kaWFnbm9zdGljcyBjYW4gYmUgY2hlY2tlZCBlbiBtYXNzZSB1c2luZyBgcGVyZm9ybWFuY2U6OmNoZWNrX21vZGVsKG1vZGVsKWAuCgotIHBvc3RlcmlvciBwcmVkaWN0aXZlIGNoZWNrIChiYXNpY2FsbHkgdXNlcyB0aGUgbW9kZWwgdG8gZ2VuZXJhdGUgbmV3IGZha2UgZGF0YSBhbmQgY29tcGFyZXMgdGhhdCBkYXRhIHRvIHRoZSBvYnNlcnZlZCBkYXRhKQotIGFzc3VtcHRpb246IGxpbmVhcml0eQotIGFzc3VtcHRpb246IGhvbW9nZW5laXR5IG9mIHJlc2lkdWFsIHZhcmlhbmNlIGFrYSBob21vc2NlZAotIGRpYWdub3N0aWNzOiBpbmZsdWVudGlhbCBvYnNlcnZhdGlvbnMgCi0gZGlhZ25vc3RpY3M6IGNvbGxpbmVhcml0eQotIGFzc3VtcHRpb246IG5vcm1hbGl0eSBvZiByZXNpZHVhbHMKCmBgYHtyfQpjaGVja19tb2RlbChtX2dvb2QpCmBgYAoKYGBge3J9CmNoZWNrX21vZGVsKG1fYmV0YSkKYGBgCgpgYGB7cn0KY2hlY2tfbW9kZWwobV9jb3JyKQpgYGAKCgpgYGB7cn0KY2hlY2tfbW9kZWwobV9pbmZsKQpgYGAKCmBgYHtyfQpjaGVja19tb2RlbChtX2xldmcpCmBgYAoKYGBge3J9CmNoZWNrX21vZGVsKG1fb3V0bCkKYGBgCgpgYGB7cn0KY2hlY2tfbW9kZWwobV91bmNvcnIpCmBgYAoKCmBgYHtyfQpjaGVja19tb2RlbChtX3VuZXF1YWwpCmBgYA==